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We study the dynamics of the Olami-Feder-Christensen (OFC) model of earthquakes, focusing on the be- 
havior of sequences of epicenters regarded as a growing complex network. Besides making a detailed and 
quantitative study of the effects of the borders (the occurrence of epicenters is dominated by a strong border 
effect which does not scale with system size), we examine the degree distribution and the degree correlation 
of the graph. We detect sharp differences between the conservative and nonconservative regimes of the model. 
Removing border effects, the conservative regime exhibits a Poisson-like degree statistics and is uncorrelated, 
while the nonconservative has a broad power-law-like distribution of degrees (if the smallest events are ignored), 
which reproduces the observed behavior of real earthquakes. In this regime the graph has also a unusually strong 
degree correlation among the vertices with higher degree, which is the result of the existence of temporary attrac- 
tors for the dynamics: as the system evolves, the epicenters concentrate increasingly on fewer sites, exhibiting 
strong synchronization, but eventually spread again over the lattice after a series of sufficiently large earth- 
quakes. We propose an analytical description of the dynamics of this growing network, considering a Markov 
process network with hidden variables, which is able to account for the mentioned properties. 
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I. INTRODUCTION 

Several different phenomena in nature spontaneously ex- 
hibit scale invariant statistics. An attempt to identify a sup- 
posed basic mechanism behind this behavior was made by 
Bak et al [1], who introduced the concept of Self-Organized 
Criticality (SOC). SOC is characterized by slowly driven sys- 
tems, with fast avalanche-like bursts of dissipation. Despite 
probably not being the sole explanation for scale-invariance 
in nature, a wide range of systems do appear to exhibit SOC, 
such as sand piles [1], forest fires [2] and earthquakes [3]. 
However, no general framework for SOC systems exist, and 
the mechanism behind it is nor very well understood. In par- 
ticular, the existence of SOC in nonconservative systems is 
still debated [4—6]. This discussion is frequently focused on 
one of the most studied and archetypal nonconservative SOC 
models, the Olami-Feder-Christensen (OFC) model for earth- 
quakes. Despite being defined by very simple rules (see sec- 
tion II), this model possesses very rich dynamics, and is able 
to reproduce a wide range of statistics of real earthquakes, 
such as the Gutenberg-Richter law for the distribution of event 
sizes [3, 7] and the Omori law for fore- and aftershocks [8, 9]. 

In this work we concentrate on the behavior of the epicen- 
ters in the OFC model, both in the conservative and noncon- 
servative regime, studied as a growing complex network with 
scale free behavior [10, 11]. 

As known previously [4, 12], we confirm that in both 
regimes epicenters are more frequent closer to the border, and 
study this effect in detail. We show, however, that this bor- 
der effect does not scale with system size, and should not 
therefore be considered representative of the dynamics of the 
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model in the thermodynamic limit. The length of the ef- 
fect is dependent on the level of dissipation, and is relatively 
large for the range of parameters normally studied, specially 
when close to the conservative limit, where a exponentially- 
decaying layer dominates, and it is hard to observe anything 
else other than this border effect. The existence of this non- 
scaling border is in accordance with what was found in [4], 
that only earthquakes from a smaller internal subset of the lat- 
tice exhibits Finite-Size Scaling in the event size statistics. 

We turn then to the dynamics of epicenters. Recently there 
has been an increasing interest in Complex Networks [13] as 
a tool for describing very diverse systems, many of which ex- 
hibit a type of scale invariance, that seems to be due to a gen- 
eral mechanism of preferential attachment [14, 15]. Inorderto 
study the epicenter dynamics in the OFC model, we construct 
a network of consecutive epicenters in the bulk, and examine 
its properties in more detail [27] . 

The network of epicenters, in the nonconservative regime, 
shows scale invariance in the degree statistics, if the epicen- 
ters of the smaller events are discarded. This network has also 
an unusual correlation among vertices of high degree, which 
makes it very distinct from networks created with a preferen- 
tial attachment rule. These results reproduce what has been 
found by Abe and Suzuki [16, 17] for real earthquakes, fur- 
ther contributing to the success of this simple model in cap- 
turing the essential earthquake dynamics. We show that this 
degree correlation seems to be due to the existence of tempo- 
rary attractors for the dynamics, that shows periods of strong 
synchronization. We also noticed that a drop in the average 
in-degree of the network seems to precede big earthquakes, 
what could in principle be used to predict at least an increase 
in the probability of big events in a given fault. 

We also show that is possible to reproduce some of the char- 
acteristics of the complex epicenter network found in the non- 
conservative version of the OFC model defining a growing 
procedure based on a Markov Chain with hidden variables. To 
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each possible epicenter (vertex) is attached a hidden variable, 
and the probability of connections among epicenters (related 
to the time sequence of events) is now given as a function 
of the hidden variable of both vertices (instead of a simple 
preferential attachment rule, as in a Barabasi- Albert type net- 
work [15]). 

This paper is organized as follows: in section II we briefly 
present the Olami-Feder-Christensen model for earthquakes; 
in section III we discuss in detail the way the spatial distribu- 
tion of epicenters depends on the distance to the borders, in 
the conservative, nonconservative and "almost conservative" 
regimes; in section IV we review the way we can built a scale 
free network from the time series of epicenters, and present 
the main properties of this network when the border effect is 
discarded. This network, although showing a scale-free be- 
havior, is quite different from Barabasi-Albert type networks, 
with a strong correlation among vertices with high degree. In 
section V we show how we can grow a network with simi- 
lar properties based on a Markov Chain process with hidden 
variables and finally, in section VI, we summarize our results. 
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II. THE OFC MODEL 

The OFC model [3] was inspired by the Burridge-Knopoff 
spring -block model [18], and is defined as a 2D coupled map 
on a square lattice. To each site in the lattice is assigned 
a "tension" Zij, initially chosen at random from the interval 
[0, z c [. The entire system is driven slowly, with every in- 
creasing uniformly. Whenever a site reaches the threshold ten- 
sion (z^ = z c ), an avalanche starts (the "earthquake"). The 
first site to reach z c and start an avalanche is called the epicen- 
ter. A fraction a of the tension of the toppling site is transfered 
to each of its four neighbors (2j±i,j±i = z i±i,j±i + otZij), 
and its tension is set to zero (z^ = 0). If any of the neigh- 
bors acquires a tension Zjiijii > z c , the same toppling rules 
are applied, until there are no more sites in the system with 
z^ > z c . Without loss of generality, we set z c = 1, The total 
number of sites that topple until the avalanche is over is called 
the "size" of the avalanche. The parameter a defines the level 
of local conservation of the system. For a = 0.25 the system 
is locally conservative and for a < 0.25 it is dissipative. We 
consider here only the case with open boundary conditions, 
i.e., the sites at the border of the lattice transfer tension to 
nonexisting neighbors, so the system is always globally non- 
conservative, but tends to conservative in the thermodynamic 
limit if a = 0.25. 



III. INFLUENCE OF THE BORDERS IN THE 
FREQUENCY OF EPICENTERS 

We find that, in the stationary regime of the OFC model, 
the number of times a site is an epicenter varies according to 
how close that site is from the border, with epicenters closer 
to the border occurring much more often. We will refer to 
this excess of epicenters in the borders as the border effect. 
Fig. 1 shows the average frequency in which a site was an 



FIG. 1: (color online) Frequency of epicenters (y axis) as a 
function of the distance from the border (x axis), for different 
values of a, L and earthquakes sizes s. The data for L = 400 
in (a) was shifted upwards for clarity. From (c) to (f) the data 
for different values of L were collapsed on top of each other 
by hand. All quantities are dimensionless. 



epicenter, given its distance from the border, for a = 0.25, 
a = 0.249, a = 0.22 and a = 0.18. We have gathered 
statistics from two lattice sizes, L = 400 and L = 800, and 
considered at least 6 x 10 6 events (after the transient). We have 
considered epicenters only from earthquakes larger than one 
(s > 2), since size one earthquakes seem to obey their own 
statistics [19]. We have also considered epicenters that gave 
rise to larger earthquakes (s > 30), to observe the dependence 
of the border effect to earthquake size [28]. 

For the conservative regime, as can be seen in Fig. 1(a), 
the border effect is clearly weaker than in the nonconservative 
regime, and the decay proceeds slowly towards the bulk. It 
does not seem to scale with system size. Moreover, the de- 
pendence on earthquake size appears to be weak for most of 
the border effect, except for the very first few layers of sites 
close to the border. 

Figures 1(c) to 1(f) also shows the same results for the non- 
conservative case, for a = 0.18 and a = 0.22. We notice 
that the border effect is composed roughly of three parts: A 
thin region, comprised of the first few sites closest to the bor- 
der, where the effect is strong and seems to decay exponen- 
tially. This region is followed by a thicker layer of sites, with 
a slower but also exponentially-decaying effect, and finally 
there is a third region in which the decay is not exponential 
and proceeds still more slowly towards the bulk of the system. 
None of this regions seems to scale with system size, with the 
possible exception of the third longer layer. The overall border 
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effect seems, however, to depend on the earthquake size (on 
the contrary of what was observed in the conservative case), as 
can be seen in Figs. 1(d) and 1(f), which shows clearly that the 
border effect decays more slowly towards the bulk of the sys- 
tem if only larger events are considered. In Figs. 1(c) to 1(f), 
the data for lattices of different size L was collapsed by hand, 
that is, curves were shifted up and down in order to coincide, 
since statistics are different in each case. The slope and the 
size of the layers, however, were not changed. 

The border effect also depends on a. The closer the system 
is to the conservative regime, the stronger and thicker is the 
layer of sites affected by it. Note that for a = 0.249, the 
"almost conservative" case (see Fig. 1(b)), the border effect 
is so strong that almost no epicenters happen in the bulk of 
the system, and only the fast exponentially-decaying border 
effect is seen. This indicates that the lattice size L = 800 
is still too small to study the system in this regime. If we 
compare this figure with Fig. 1(a), we note that there is also 
an evidence of a sharp transition from the nonconservative to 
the conservative regimes of the model, for which the border 
dependence of epicenters is radically different. 

The crucial role of the border in this model was already 
pointed out by Middleton and Tang [12], who argued that the 
inhomogeneity introduced by the open boundary inhibits the 
synchronization of the bulk, which would otherwise reach a 
periodic state, as it happens with the system with periodic 
boundaries. The resulting "self-organization" would begin 
at the border and then proceed towards the bulk, following a 
power-law in time. Also, it has been shown in [4], that while 
the statics of event sizes in the OFC model does not seem to 
obey Finite Size Scaling (FSS), this behavior is recovered only 
when events inside a smaller internal subset of the tension lat- 
tice are considered. Thus, the existence of non-scaling border 
effects is to be expected. 

We proceed to examine the dynamics of the epicenters in 
the model, but only those unrelated to the non-scaling border 
effect. Therefore, unless otherwise noted, we ignored all the 
epicenters belonging to an outer layer of 100 sites in the lat- 
tice, for all the systems studied. 



IV. SEQUENCES OF EPICENTERS AS A COMPLEX 
NETWORK 

A graph (or network) is a set of discrete items, called ver- 
tices or nodes, with connections between them, called edges 
or links. An edge, connecting vertices i and j, is directed if it 
is defined in only one direction (connects vertex i with vertex 
j, for instance, but not site j with site i) and a graph is said 
to be directed if its edges are directed. There may be more 
than one edge between a pair of vertices, and the graph is 
called in this case a multigraph. The number of edges con- 
nected to a vertex is called the degree of the vertex; since 
there may be more than one edge between two vertices, the 
degree of a vertex is not necessarily equal to the number of 
its neighbors. If the graph is directed, it is then possible to 
talk about out-degree (number of edges leaving a vertex) and 
in-degree (number of edges incident to a vertex). The degree 




FIG. 2: (color online) An example of an epicenter graph 
(right), generated from a sequence of epicenters (marked in 
red) from a 5x5 tension lattice (left). The graph corresponds 
to the following sequence of epicenters: 3, 20, 23, 13, 20, 17, 
2. 

distribution of a graph, P(k), gives the probability that a ran- 
domly sampled vertex has degree k. Graphs have been the 
subject of systematic study by mathematicians for some time, 
but recent years witnessed a growth in the interest on this sub- 
ject among physicists, with emphasis on large-scale statistical 
properties of graphs. Many statistical mechanics concepts and 
techniques have been widely used, and a good review on re- 
cent developments in this subject can be found in [13]. We 
will show that some tools of network theory can help to get a 
deeper understanding of the dynamics of the OFC model and 
maybe of the dynamics of real earthquakes. 

The sequence of epicenters in the OFC model can be used to 
construct a directed multigraph in the following manner. Each 
site that is an epicenter represents a vertex. Two consecutive 
epicenters are connected by a directed edge, from the first to 
occur to the second (see Fig. 2). Since, in principle, the same 
site can become an epicenter two times consecutively, loops 
are allowed (but don't occur often). It is also possible for the 
same sequence of epicenters to happen more than once, so 
parallel edges are also allowed. This graph has certain reg- 
ularities: The out-degree of every vertex is always equal to 
the in-degree, except for the very first and last epicenters of 
the sequence, and therefore the total degree is always an even 
number. Also, if the direction of the edges is ignored, the 
graph is always composed of only one component. 

We have constructed graphs for the epicenters of the OFC 
model with L = 400 and 800, and for a = 0.25, a = 0.22 
and a = 0.18. We also considered the graphs for epicenters 
of different earthquakes sizes. We then observed the degree 
distribution and the degree correlation of the graph. The re- 
sults for the nonconservative regime are averages over 5 to 1 1 
graphs, depending on the size of earthquakes considered, each 
graph with 6 x 10 6 edges. 



A. Degree distribution 

Since the in-degree of the network is equal to the out- 
degree, it is sufficient to describe only one of the two, and 
here we choose arbitrarily the in-degree. 

For the conservative regime (Fig. 3(a)), the in-degree dis- 
tribution seems to be a Poisson (which gets stretched if more 
sites from the border are considered), indicating that, in this 
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FIG. 3: (color online) In-degree distribution P(k) (y axis) in 
function of the in-degree k (x axis), for different values of L, 
a, and earthquake sizes s. In (a) are shown the distributions 
for two different sizes of the discarded border b, and the solid 
line is the corresponding Poisson distribution. The data for 
b = was shifted upwards for clarity. In (b) and (c) the solid 
line is the result of fitting a power-law k 'to the data when 
s > 30. The data for different earthquake sizes were shifted 
upwards for clarity. All quantities are dimensionless. 



regime, epicenters in the bulk of the lattice occur randomly. 
Moreover, the degree distribution does not depend on the min- 
imum size of the earthquakes considered. 

For the nonconservative regime the situation changes. As 
can be seen in Fig. 3(b) and 3(c), if only larger earthquakes 
are considered, the in-degree distribution resembles more a 
power-law. The exponent of the power-law seem to be de- 
pendent on a, with smaller a leading to steeper lines. For 
a = 0.22 and s > 30, in Fig. 3(c), the high fluctuations at 
the tail of the in-degree distribution represent a lack of statis- 
tics, due to an average over only five realizations of the graph, 
while for s > 5, for instance, the average was over ten differ- 
ent graphs. For both a = 0.22 and a = 0.18, the difference 
of inclination of the power-law region of the distributions is 
very small between the data for s > 5 and s > 30, indicat- 
ing that it is not strongly dependent on the lower bound of the 
considered earthquake sizes, provided it is large enough for 
the power-law to emerge. 



Correlations between degree distribution and tension 
distribution in the lattice 



It is interesting to observe where the epicenters happen in 
the tension lattice. As has already been shown in [20], the 
stationary state of the OFC model, for a < 0.25 (nonconser- 
vative), exhibits patchy synchronized regions within the bulk 
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FIG. 4: (color online) Snapshot of the tension lattice at the 
stationary state, for L = 800 and a = 0.18. The next 10 4 
epicenters, for earthquake sizes s > 2, after this configuration, 
are marked in green. All quantities are dimensionless. 



of the system with sites that have similar tension, and behave 
similarly to the OFC model with periodic boundary condi- 
tions, exhibiting heavy synchronization. As can be seen in 
Fig. 4, for a = 0.18, the epicenters seem to happen mostly in 
the frontiers among those synchronized regions, and in valley- 
like structures inside the plateaus. As only larger earthquakes 
are considered, the epicenters happen increasingly in smaller 
and less structured regions (not shown). The same behavior 
was also observed for a = 0.22. 

In Fig. 5 can be seen the in-degree of a vertex placed in 
the tension lattice, i.e., the number of times a site was an epi- 
center, for a = 0.18 and s > 2. The epicenters seem to be 
well distributed inside the bulk, but aggregated in stripe-like 
structures. For s > 30 the epicenters seem considerably less 
aggregated (not shown). For a = 0.22 the results were ob- 
served to be very similar. 



C. Degree correlation 

One further basic aspect of the epicenter network which we 
analyzed was the degree correlation, i.e., how vertices are con- 
nected to each other based on their degrees. We look at the 
average in-degree of the nearest "out-neighbors" of a vertex 
(vertices that receive an edge coming from it), in function of 
the degree of the vertex. 

We found that for the conservative regime (Fig. 6(a)), the 
graph seems to be uncorrelated, with the in-degree of the near- 
est neighbors being independent on the in-degree of the origi- 
nating vertex. Together with the in-degree distribution (a Pois- 
son), this puts this graph closer to the class of totally random 
graphs such as the Erdos-Renyi graph [21]. 

The situation is again very different for the nonconservative 
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FIG. 5: (color online) In-degree of vertices placed in the bulk 
of tension lattice, i.e., the number of times a site was an epi- 
center, for L = 800 and a = 0.18. Only earthquakes with 
sizes s > 2 were considered. All quantities are dimension- 
less. 
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FIG. 6: (color online) Average in-degree of nearest out- 
neighbors k nn (k) (y axis) in function of the in-degree k (x 
axis), for different values of L, a, and earthquake sizes s. The 
solid lines are fitted straight lines. All quantities are dimen- 
sionless. 



regime (Figs. 6(b) and 6(c)). In that case, the degrees seem 
highly correlated, with vertices with high in-degree connect- 
ing predominantly to other vertices of high in-degree, which 
makes the network assortative. The correlation seems to be 
linear for higher degrees, when only larger earthquakes are 
considered [29]. Citation networks [14, 15] and other net- 
works that are grown with a preferential attachment rule have 



a quite different behavior, with an in-degree distribution fol- 
lowing a power-law, but in those cases the degree correlation 
also decays with a power law [22], converging to a constant 
value for large in-degrees. Thus, the dynamics responsible for 
generating this network must be fundamentally different than 
the dynamics generated by a preferential attachment rule. Re- 
cently it has also been found that a very similar network, when 
constructed with real earthquake data, is also assortative and 
exhibits similar degree correlations [17]. 

What indeed is unveiled by this high correlation amongst 
high in-degree vertices is an attracting dynamics: Connec- 
tions from vertices of one type are much more probable to 
vertices of the same type, eventually trapping the sequence 
of epicenters in a smaller region of the lattice, stretching the 
in-degree distribution, and generating the observed in-degree 
correlation. This trapping seems to be strongly correlated to 
the occurrence of very large earthquakes, and the large scale 
redistribution of tensions that is caused by them. This can be 
seen in Fig. 7, where is shown the average in-degree of the 
subgraph composed only of the last 10 5 events, together with 
the amplitude of the corresponding events. Whenever a large 
earthquake occurs, the average in-degree drops, meaning that 
the last epicenters happened in a larger number of sites. In 
fact the decay of the average in-degree starts before the main 
big earthquake, and seems to occur together with the smaller 
events that lead up to it, the so called fore shocks [8, 9, 23]. 
Thus, the large events, together with their foreshocks, are re- 
sponsible for breaking the attractor, and spreading the epicen- 
ters to a larger region. After the sequence of large events, the 
trapping of epicenters starts again, until the next sequence of 
large events sweeps it again. Although we did not make an ex- 
tensive analysis to define the degree of certitude of the this ob- 
servation, monitoring the in-degree of this network may repre- 
sent a promising way of predicting an increase in the probabil- 
ity of observing large earthquakes in a given fault, and to iden- 
tify, among the small events, the signature of the foreshocks 
that preceed a main shock [30]. Since the network of epicen- 
ters generated by the OFC model seems to reproduce many as- 
pects of the network of epicenters built from real data [16, 17], 
including the degree correlation mentioned above, it would be 
interesting to see in more detail if both graphs are actually 
generated by the same overall dynamics. This however would 
require a more systematic and thorough analysis of real earth- 
quake data, and therefore would be better suited for a separate 
work. 

To illustrate the topology of the graph during both situa- 
tions, we show a subgraph of the whole network, correspond- 
ing to a region of 10 4 events collected during the period that 
the dynamics is trapped in an attractor (Fig. 8), and just af- 
ter a large earthquake (Fig. 9), as indicated in Fig. 7. As can 
be seen in Fig. 8, the attractor region is dominated by syn- 
chronization, where the same sequence of ~ 10 epicenters 
occur repeatedly. During the occurrence of the large events, 
the same subgraph looks like Fig. 9, where synchronization is 
still present, but in a much smaller degree. 
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FIG. 7: (color online) Average in-degree of the subgraph com- 
posed only of the last 10 5 events, and the amplitudes of the 
events that generated the graph. The regions indicated by the 
arrows correspond to the subgraphs in Figs. 8 and 9. All quan- 
tities are dimensionless. 




FIG. 8: (color online) Subgraph composed of 10 4 consecutive 
epicenters, corresponding to the marked region at the right in 
Fig. 7. 



V. MARKOV NETWORKS WITH HIDDEN VARIABLES 

In this section we describe a general random graph model, 
based on hidden variables and a Markov Chain. It is based on 
a similar class of networks developed by Bogufia et al [24], but 
modified in order to account for the topology of the epicenter 
graph observed in the OFC model. Our goal is to better un- 
derstand the type of dynamics that is able to generate graphs 




FIG. 9: (color online) Subgraph composed of 10 4 consecutive 
epicenters, corresponding to the marked region at the left in 
Fig. 7. 

with properties examined in the previous section. 

Consider a set of N vertices, where N S> 1. To each vertex 
v is assigned a hidden variable h v , sampled from a distribution 
p(h). A directed multigraph can be constructed by a Markov 
chain, in the following manner: Starting from a random ver- 
tex /a, a directed edge is added from ptov with probability 
P(fj, —>i>)= rihfifhv), and likewise from v to any other 
vertex co with probability given by r{h v , h u ), and so forth. 
After a transient stage, the graph will have properties that are 
entirely defined by p(h) and r(h, h'). This graph is rather 
general, and, in fact, every Markov process generates such a 
graph if the discrete states of the chain are thought as vertices 
and the transition as directed edges. With this basic procedure 
in mind we can proceed to calculate the statistical properties 
of the graph. 

A. In-degree distribution 

Like in the network of epicenters, every vertex of the net- 
work generated in the way described above has the in-degree 
equal to the out-degree. Thus, it is sufficient to describe only 
one of the two. To find the in-degree distribution of this graph, 
one must consider an ensemble of graphs and the probability 
in the ensemble of one vertex v receiving one connections af- 
ter a time T, w v (T), which is given by 

w v (T + l)='£,P(ji-n/)w lll (T). (1) 

After a long time T, the system reaches the stationary state, 
w(oo), given by 

iy(oo) = P"w(oc), (2) 

where P is the transition matrix defined by P(fx — > v), w(T) 
is the state vector at time T, and n is the period of the solution 
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(we will consider only n = 1 from now on). 

The probability that a vertex p has in-degree fc after a time 
T> 1, P(k\p, T), is given simply by the binomial distribu- 
tion, 



k„-Tw u 



where Wk = Wfc (oo), which can be approximated by the Pois- 
son distribution, as in the rightmost term. 

The total in-degree distribution after a time T, P(k\T), is 
then given by 



(4) 



Now since a vertex p is labeled uniquely by its hidden vari- 
able h,j,, we must have then that = w(h u ). Thus, w(h) 
can be obtained by rewriting equation 1, 



w(h) = N J r(h tl ,h)w(h f _ t )p(h t _ t )dh f _ t , (5) 

h 

assuming that h is a continuous variable (the last expression 
would just be a sum if it were discrete). Solving this inte- 
gral equation for w(h), it is possible then to obtain the degree 
distribution through equation 4, 



1. In-degree correlation 

It is also possible to calculate the degree correlation of this 
graph. The probability of one vertex p, with in-degree k, con- 
necting to another vertex of degree k' is given by 

P(fc ' |fc ' ^ T) = P NP{k T T) ^ Pijl " V)P[k ' ~ lk ' T) - 

(7) 

The total probability of one vertex with degree k connecting 
to one of degree k' is then simply 



P{k'\k,T)=Y,P{k'\k, l i,T), 



(8) 



and the average in-degree of the nearest out-neighbors is just 
then 



h m {k,T) = Y,k'P{k'\k,T). 



(9) 



In terms of the hidden variables, substituting equations 6 and 3 
in 8 and calculating the sum in 9, we have then, 



fcnn(fc,T) = 1 + 



N f f (Tw(h^)) k e - T ™(M 



kl 



P(k,T) 

h 

x r(h l _ l ,K)Tw(h v )p(h l _ l )p(h^)dh fJ ,dh u . (10) 









" hi 


G{h v ) 1 




p(. h fj.) J 






hv 




1 >- 



h u ) x p(h v ) 



Kp{h v ) 




(a) 



FIG. 10: (a) Connection probability (equation 11) from a ver- 
tex p, to a vertex v, and (b) connection probability from a ver- 
tex p, to any vertex with hidden variable h v . 



p(h) 


G{h) 




P(k) 


h-f 
p-1 


1 


~ k 


~ k 7+i 




1 


~ k 








~ k 


~ fc ~ 



TABLE I: Different asymptotic shapes for k nn (k) and P(k) 
for different shapes of G(h) and p(h), for fc 3> 1. 



2. Attractor dynamics 

We want to understand how correlations such as seen in 
Figs. 6(b) and 6(c), and power-law distributions can arise 
from this type of network. For that we must define a suitable 
r(h, h') and p(h). It is clear that what uniquely defines the 
in-degree of some vertex is its hidden variable. Thus, for the 
in-degree correlation to be of the form fc nn (fc) ~ k for large 
fc, we must have that h'(h) ~ h for large h, where h'(h) is 
the average hidden variable of the out-neighbors of a vertex 
with hidden variable h. With this in mind, we define then the 
following general expression for the connection probability, 



r(h^,h v )=F(h l ,) 



G{K)hl 
p(h v ) 



[h v < hp] 



[K > h. 



(11) 



Where G(h) is a function that dictates how fast the connec- 
tion probability decays for h v < (see Fig. 10), and the 
exponent 7 defines the preference with which vertices with 
higher h are chosen. The function F(h^,) is simply given by 
the normalization condition J2 U r (^^, h v ) = 1. 

We considered a few shapes for G(h) and p(h) and calcu- 
lated the degree distribution and degree correlation through 
equations 6 and 10, always for fc ;§> 1. The results are sum- 
marized in table I. 

What we find is that the effect of adopting a connection 
probability like the one described by equation 1 1 is to gener- 
ate a in-degree distribution corresponding to a stretched form 
of p(h). If G(h) is independent of h, and the trapping in the 
region of similar h is the weakest, we have the following pos- 
sibilities: If p(h) is a power-law with exponent j3, then P(fc) 



8 



will also be a power-law with exponent in the region [1,(3], 
approaching 1 if 7 is large. When p(h) is a exponential distri- 
bution, the resulting in-degree distribution will be a stretched 
exponential as indicated in table I, which will also resemble a 
power-law if 7 is relatively large. Now, considering a stronger 
"trapping effect" with G(h) increasing exponentially, we have 
that an exponential p(h), with decay parameter (3, is enough to 
create a power-law distribution of in-degrees, with exponents 
in the interval [l,j3], approaching 1 with faster G(h). This 
means that it is not necessary to assume an intrinsic scale- 
invariance, represented by a power-law in p(h), for the exis- 
tence of a power-law in P(k). Furthermore, the asymptotic 
in-degree distribution in this case does not depend on 7, being 
totally dominated by the "trapping" behavior, and not by the 
preference of connection. 

The process described above shows a variety of ways in 
which graphs with in-degree distributions resembling power- 
laws and linear in-degree correlation can be created. Looking 
at only these properties, it is not possible to know which one of 
the possibilities (if any) is more likely to describe the epicenter 
network. Moreover, the process above would not account for 
the strong synchronization observed in Figs. 8 and 9. After all, 
the sequence of epicenters are probably not simple Markovian 
processes. However, the above model, as a first approxima- 
tion, serves the purpose of illustrating how such correlations 
and in-degree distribution can come to place, and presents a 
general analytical framework for further modeling. 

VI. CONCLUSIONS 

We have shown that the epicenters in the OFC model oc- 
cur predominantly near the boundary of the lattice, but this 
preference does not seem to scale with system size. This 
border affinity depends on the dissipation parameter a, be- 
ing thinner for smaller values of a. It is also dependent on the 
earthquake size, with epicenters of larger earthquakes having 
a border effect which decays more slowly towards the bulk. 
We have also studied the network of consecutive epicenters, 
and found that it is sharply different in the two regimes of 
the model: In the conservative regime it is rather feature- 
less, with uncorrelated in-degree statistics and Poisson in- 



degree distribution. However, in the nonconservative regime, 
it has an unusual linear degree correlation amongst vertices 
of high degree, and a broad distribution of in-degrees resem- 
bling a power-law, but only when the smaller earthquakes are 
not considered. The in-degree distribution and correlation in 
this regime is similar to what was found very recently for 
real earthquakes [16, 17]. Furthermore, we noticed that the 
high correlation of in-degrees is due to an attractor dynamics 
where the occurrence of epicenters tend to synchronize, with 
the same sequence of epicenters occurring continuously. This 
synchronization is broken by large earthquakes, which spread 
the epicenters over a larger portion of the lattice, thus populat- 
ing the graph with vertices of smaller in-degree. Interestingly, 
the effects of the large events on the topology of the epicenter 
network are noticeable before the actual main event, and seem 
to be related to the series of increasingly larger foreshocks 
that precede it. Since the prediction of the OFC model that 
there would be an in-degree correlation in the epicenter graph 
corresponds to what has been recently found for real earth- 
quakes [17], further detailed analysis of this behavior may 
prove useful for the prediction of large earthquakes. Lastly 
we described a general analytical network model based on 
a Markovian process and hidden variables, which is able to 
reproduce the most general aspects of the epicenter network, 
when a suitable attractor dynamics is specified. There are sev- 
eral aspects of the dynamics of epicenters that remain uncov- 
ered. It would be of special interest to look at other topologi- 
cal properties of the epicenter graph, such as the dependence 
of the clustering coefficient on in-degree, and the existence of 
community structure [25, 26]. Furthermore it would also be 
useful to compare in detail some of the results here obtained, 
such as the dynamics responsible for the in-degree correlation 
and the epicenter synchronization, with the epicenter network 
of real earthquakes. 
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